## Descriptive statistics
## Brian T. Hamel and Bryan Wilcox-Archuleta
## First: 28 April 2019  
## Last: 19 March 2020

## Loading packages
## install.packages(c("tidyverse", "estimatr", "ggpubr", "gridExtra", "texreg", "scales"))
library(tidyverse)
library(estimatr)
library(ggpubr)
library(gridExtra)
library(texreg)
library(scales)

## Loading Racial Flux data
load("01_data/lodes/racial_flux.RData")

## Loading ACS data
load("01_data/acs/acs.RData")

## Merging Racial Flux and ACS data
dta = racial_flux %>%
  left_join(., acs, by = c("zcta" = "zip"))

## Plotting Racial Flux vs. % white and % black
###############
## FIGURE 1 ##
###############
pct_white = ggplot(dta, aes(x = pct_white, y = racial_flux)) + 
  geom_point(color = "black", shape = 1, alpha = .25) + 
  stat_smooth(method = "lm_robust", se = FALSE, lty = 2, color = "blue", show.legend = TRUE) +
  stat_smooth(method = "loess", se = FALSE, color = "red", show.legend = TRUE) + 
  labs(x = "% White", y = "Racial Flux") + 
  scale_y_continuous(labels = number_format(accuracy = 0.01)) +
  theme(legend.title = element_blank()) +
  stat_cor(method = "pearson", label.x = 2, label.y = 85, size = 2)

pct_black = ggplot(dta, aes(x = pct_black, y = racial_flux)) + 
  geom_point(color = "black", shape = 1, alpha = .25) + 
  stat_smooth(method = "lm_robust", se = FALSE, lty = 2, color = "blue", show.legend = TRUE) +
  stat_smooth(method = "loess", se = FALSE, color = "red", show.legend = TRUE) + 
  labs(x = "% Black", y = "") +
  scale_y_continuous(labels = number_format(accuracy = 0.01)) +
  theme(legend.title = element_blank()) +
  stat_cor(method = "pearson", label.x = 2, label.y = 85, size = 2)

## Combine plots, and save
flux_res = grid.arrange(pct_white, pct_black, ncol = 2, nrow = 1)
ggsave(flux_res, file = "03_output/flux_res.png", height = 2, width = 4, units = "in", dpi = 600)

## Correlates of Racial Flux
correlates = lm_robust(racial_flux ~ pct_white + pct_black
                       + pct_unemployed + pct_college + log_per_cap_inc + gini + south
                       + non_rural + log_pop_density, data = dta, se_type = "stata")

## Save table of coefficients
##############
## TABLE A2 ##
##############
texreg(correlates,
       file = "03_output/correlates.tex", 
       label = "correlates",
       caption = "Multivariate Correlates of Racial Flux", 
       custom.model.names = c("(1)"),
       custom.coef.names = c("Intercept", "% White", "% Black", 
                             "% Unemployed", "% College", 
                             "log(Per Capita Income)", "Gini Coef.", "South",
                             "Non-Rural", "log(Pop. Density)"), 
       reorder.coef = c(2, 3, 4, 5, 6, 7, 8, 9, 10, 1),
       custom.gof.names = c(NA, NA, "Observations", NA, NA),
       stars = c(0.05, 0.01, 0.001),
       digits = 3,
       center = TRUE,
       include.ci = FALSE,
       caption.above = TRUE)

## Clear R 
rm(list = ls())
